home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
fortran
/
linpkdrv.zip
/
ST.FOR
< prev
next >
Wrap
Text File
|
1984-01-12
|
16KB
|
546 lines
C MAIN PROGRAM
INTEGER LUNIT
C ALLOW 5000 UNDERFLOWS.
C CALL TRAPS(0,0,5001,0,0)
C
C OUTPUT UNIT NUMBER
C
LUNIT = 6
C
CALL STRTS(LUNIT)
STOP
END
SUBROUTINE STRTS(LUNIT)
INTEGER LUNIT
C LUNIT IS THE OUTPUT UNIT NUMBER
C
C TESTS
C STRCO, STRSL
C
C LINPACK. THIS VERSION DATED 08/14/78 .
C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C SUBROUTINES AND FUNCTIONS
C
C LINPACK STRCO,STRSL
C EXTERNAL STRXX,SMACH
C BLAS SAXPY,SDOT,SSCAL,SASUM
C FORTRAN ABS,AMAX1,FLOAT,MAX0
C
C INTERNAL VARIABLES
C
REAL A(15,15),B(15),BT(15),X(15),XEXACT(15),XT(15),Z(15)
REAL AINV(15,15),DET(2)
REAL SDOT,STUFF,T
REAL ANORM,AINORM,RCOND,COND,COND1,SMACH,EPS
REAL ENORM,ETNORM,RNORM,RTNORM,XNORM,XTNORM,EN,SASUM
REAL FNORM,ONEPX,Q(7),QS(7)
INTEGER I,INFO,J,JOB,KASE,KOUNT,KSING,LDA,ML,MU,N,NPRINT
INTEGER KSUSP(7),IQ(7)
C
LDA = 15
C
C WRITE MATRIX AND SOLUTIONS IF N .LE. NPRINT
C
NPRINT = 3
C
WRITE (LUNIT,380)
WRITE (LUNIT,730)
C
DO 10 I = 1, 7
KSUSP(I) = 0
10 CONTINUE
KSING = 0
C
C SET EPS TO ROUNDING UNIT
C
EPS = SMACH(1)
WRITE (LUNIT,390) EPS
WRITE (LUNIT,370)
C
C START MAIN LOOP
C
KASE = 1
20 CONTINUE
C
C GENERATE TEST MATRIX
C
CALL STRXX(A,LDA,N,KASE,LUNIT)
C
C N = 0 SIGNALS NO MORE TEST MATRICES
C
C ...EXIT
IF (N .LE. 0) GO TO 360
ANORM = 0.0E0
DO 30 J = 1, N
ANORM = AMAX1(ANORM,SASUM(N,A(1,J),1))
30 CONTINUE
WRITE (LUNIT,570) ANORM
C
IF (N .GT. NPRINT) GO TO 50
WRITE (LUNIT,370)
DO 40 I = 1, N
WRITE (LUNIT,600) (A(I,J), J = 1, N)
40 CONTINUE
WRITE (LUNIT,370)
50 CONTINUE
C
C GENERATE EXACT SOLUTION
C
XEXACT(1) = 1.0E0
IF (N .GE. 2) XEXACT(2) = 0.0E0
IF (N .LE. 2) GO TO 70
DO 60 I = 3, N
XEXACT(I) = -XEXACT(I-2)
60 CONTINUE
70 CONTINUE
C
C GENERATE R.H.S.
C
DO 90 I = 1, N
B(I) = 0.0E0
BT(I) = 0.0E0
DO 80 J = 1, N
B(I) = B(I) + A(I,J)*XEXACT(J)
BT(I) = BT(I) + A(J,I)*XEXACT(J)
80 CONTINUE
X(I) = B(I)
XT(I) = BT(I)
90 CONTINUE
C
C UPPER OR LOWER TRIANGULAR
C
ML = 0
MU = 0
DO 120 J = 1, N
DO 110 I = 1, N
IF (A(I,J) .EQ. 0.0E0) GO TO 100
IF (I .LT. J) MU = MAX0(MU,J-I)
IF (I .GT. J) ML = MAX0(ML,I-J)
100 CONTINUE
110 CONTINUE
120 CONTINUE
WRITE (LUNIT,670) ML,MU
IF (ML .NE. 0 .AND. MU .NE. 0) GO TO 350
IF (MU .EQ. 0) JOB = 0
IF (ML .EQ. 0) JOB = 1
IF (JOB .EQ. 0) WRITE (LUNIT,710)
IF (JOB .EQ. 1) WRITE (LUNIT,720)
STUFF = 4095.0E0
DO 140 J = 1, N
DO 130 I = 1, N
IF (I .LT. J .AND. JOB .EQ. 0) A(I,J) = STUFF
IF (I .GT. J .AND. JOB .EQ. 1) A(I,J) = STUFF
130 CONTINUE
140 CONTINUE
C
C ESTIMATE CONDITION
C
CALL STRCO(A,LDA,N,RCOND,Z,JOB)
C
C OUTPUT NULL VECTOR IF N .LE. NPRINT
C
IF (N .GT. NPRINT) GO TO 160
WRITE (LUNIT,620)
DO 150 I = 1, N
WRITE (LUNIT,630) Z(I)
150 CONTINUE
WRITE (LUNIT,370)
160 CONTINUE
C
C
C TEST FOR SINGULARITY
C
IF (RCOND .GT. 0.0E0) GO TO 170
WRITE (LUNIT,610) RCOND
WRITE (LUNIT,400)
KSING = KSING + 1
GO TO 340
170 CONTINUE
COND = 1.0E0/RCOND
WRITE (LUNIT,420) COND
ONEPX = 1.0E0 + RCOND
IF (ONEPX .EQ. 1.0E0) WRITE (LUNIT,410)
C
C COMPUTE INVERSE, DETERMINANT AND COND1 = TRUE CONDITION
C
DO 190 J = 1, N
DO 180 I = 1, N
AINV(I,J) = A(I,J)
180 CONTINUE
190 CONTINUE
CALL STRDI(AINV,LDA,N,DET,110+JOB,INFO)
AINORM = 0.0E0
DO 200 J = 1, N
IF (JOB .EQ. 0)
* AINORM = AMAX1(AINORM,SASUM(N-J+1,AINV(J,J),1))
IF (JOB .EQ. 1)
* AINORM = AMAX1(AINORM,SASUM(J,AINV(1,J),1))
200 CONTINUE
COND1 = ANORM*AINORM
WRITE (LUNIT,430) COND1
WRITE (LUNIT,650) DET(1)
WRITE (LUNIT,660) DET(2)
C
C SOLVE A*X = B AND TRANS(A)*XT = BT
C
CALL STRSL(A,LDA,N,X,JOB,INFO)
CALL STRSL(A,LDA,N,XT,JOB+10,INFO)
C
IF (N .GT. NPRINT) GO TO 230
WRITE (LUNIT,440)
DO 210 I = 1, N
WRITE (LUNIT,640) X(I)
210 CONTINUE
WRITE (LUNIT,450)
DO 220 I = 1, N
WRITE (LUNIT,640) XT(I)
220 CONTINUE
WRITE (LUNIT,370)
230 CONTINUE
C
C RESTORE ZEROS IN OTHER TRIANGLE
C
DO 260 J = 1, N
DO 250 I = 1, N
IF (A(I,J) .NE. STUFF) GO TO 240
A(I,J) = 0.0E0
AINV(I,J) = 0.0E0
240 CONTINUE
250 CONTINUE
260 CONTINUE
C
C COMPUTE ERRORS AND RESIDUALS
C E = X - XEXACT
C ET = XT - XEXACT
C R = B - A*X
C RT = BT - A*XT
C AI = A*INV(A) - I
C
XNORM = SASUM(N,X,1)
XTNORM = SASUM(N,XT,1)
ENORM = 0.0E0
ETNORM = 0.0E0
FNORM = 0.0E0
DO 270 J = 1, N
ENORM = ENORM + ABS(X(J)-XEXACT(J))
ETNORM = ETNORM + ABS(XT(J)-XEXACT(J))
T = -X(J)
CALL SAXPY(N,T,A(1,J),1,B,1)
BT(J) = BT(J) - SDOT(N,A(1,J),1,XT,1)
270 CONTINUE
RNORM = SASUM(N,B,1)
RTNORM = SASUM(N,BT,1)
C
C
C A*INV(A) - I
C
AINORM = 0.0E0
DO 300 J = 1, N
DO 280 I = 1, N
B(I) = 0.0E0
280 CONTINUE
DO 290 K = 1, N
T = AINV(K,J)
CALL SAXPY(N,T,A(1,K),1,B,1)
290 CONTINUE
B(J) = B(J) - 1.0E0
AINORM = AMAX1(AINORM,SASUM(N,B,1))
300 CONTINUE
C
WRITE (LUNIT,460) ENORM,ETNORM
WRITE (LUNIT,470) RNORM,RTNORM
WRITE (LUNIT,580) AINORM
C
C COMPUTE TEST RATIOS
C
Q(1) = COND/COND1
Q(2) = COND1/COND
Q(3) = ENORM/(EPS*COND*XNORM)
Q(4) = ETNORM/(EPS*COND*XTNORM)
Q(5) = RNORM/(EPS*ANORM*XNORM)
Q(6) = RTNORM/(EPS*ANORM*XTNORM)
Q(7) = AINORM/(EPS*COND)
WRITE (LUNIT,370)
WRITE (LUNIT,480)
WRITE (LUNIT,370)
WRITE (LUNIT,540)
WRITE (LUNIT,550)
WRITE (LUNIT,560)
WRITE (LUNIT,370)
WRITE (LUNIT,590) (Q(I), I = 1, 7)
WRITE (LUNIT,370)
C
C LOOK FOR SUSPICIOUS RATIOS
C
QS(1) = 1.0E0 + 4.0E0*EPS
QS(2) = 10.0E0
EN = FLOAT(N)
IF (N .EQ. 1) EN = 2.0E0
DO 310 I = 3, 7
QS(I) = EN
310 CONTINUE
KOUNT